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1.  Introduction 


During  the  last  three  decades  meteorological  researchers  at  the  U.S.  Army  Research  Laboratory 
(ARL),  formerly  the  Atmospheric  Sciences  Laboratory  (ASL),  have  conceived,  developed,  and 
used  a  high  resolution  wind  (HRW)  model.  The  HRW  model  is  a  micro-alpha  scale  (Orlanski, 
1975),  two-dimensional,  diagnostic  model  that  simulates  wind  and  temperature  fields  in  the 
atmospheric  surface  layer,  taking  into  account  both  complex  terrain  topography  and  thermal 
structure  over  a  limited  area.  The  computational  domain  size  can  range  from  a  2  km  square  to  a 
20  km  square  with  grid  resolutions  varying  from  40  m  to  400  m,  respectively.  The  vertical 
thickness  of  the  computational  layer  is  designed  to  be  one  tenth  the  magnitude  of  the  grid  size. 

A  typical  grid  size  of  100  m,  therefore,  produces  simulated  fields  at  a  level  10  m  above  the 
surface. 

The  HRW  model  was  originally  fonnulated  by  Ball  and  Johnson  (1978)  who  described  its 
theoretical  basis  and  computational  structure  in  a  technical  report  submitted  to  ASL.  This 
geographically  re-locatable  model  was  designed  to  be  incorporated  in  the  U.S.  Army 
Experimental  Prototype  Automatic  Meteorological  System  for  the  estimation  of  the  surface  layer 
wind  field  at  sub-mesoscale  resolution  over  a  limited  area  in  broken  topography.  Later  the  HRW 
model  was  further  developed  as  a  stand-alone  model,  as  well  as  an  integral  component  of  a 
hierarchy  of  nested  meso-  and  micro-meteorological  models  (Cionco,  1987).  A  distinctive 
feature  of  the  HRW  model  is  that  it  has  adopted  a  warped  coordinate  system  to  address  the 
intimate  interaction  between  the  surface  layer  and  the  variable  ground  features.  As  far  as  the 
authors  are  aware,  this  approach  is  unique;  virtually  no  meso-  or  micro-scale  models  have  been 
formulated  with  this  approach.  The  HRW  model  has  been  tested  and  used  for  a  variety  of 
applications,  e.g.,  Weber  et  al.  (1995),  Thykier-Nielsen  et  al.  (1995),  Cionco  and  Byers  (1995, 
1997),  Cionco  (1998),  Cionco  et  al.  (1998),  Cionco  and  Ellefsen  (1998).  Although  the  model  has 
been  tested  and  used  for  a  long  time,  it  has  not  been  thoroughly  evaluated  due  to  the  fact  that 
adequate  observational  data  at  a  resolution  of  100  m  or  so  are  extremely  scarce.  Fortunately,  the 
multinational,  high-resolution  field  study  of  Meteorology  and  Diffusion  over  Non-Uniform 
Areas  (MADONA)  during  September  and  October  1992  has  provided  a  valuable  observational 
dataset.  Cionco  et  al.  (1999)  have  given  a  comprehensive  description  of  the  MADONA  project, 
including  the  dataset.  The  MADONA  field  study  was  designed  and  conducted  for,  among  other 
purposes,  high-resolution  meteorological  data  collection  in  an  effort  to  obtain  terrain-influenced 
meteorological  data.  Thirty-one  days  of  meteorological  data  were  collected.  High-resolution 
and  standard  micro-scale,  boundary  layer,  and  synoptic  meteorological  sensors  including  15 
wind  speed/direction  sets  were  deployed  over  the  MADONA  topography  (a  9  km  by  7.5  km 
area).  This  well-documented  database  is  suitable  and  valuable  for  the  evaluation  and  validation 
of  the  HRW  model. 
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The  objectives  of  this  report  are  twofold.  The  first  objective  is  to  provide  a  concise  description 
of  the  HRW  model  as  systematically  and  completely  as  possible.  The  lengthy  paper  by  Ball  and 
Johnson  (1978)  was  not  published  in  the  open  literature  and  is  not  readily  available  for  interested 
readers.  A  limited  number  of  papers  presented  in  a  book  (Cionco,  1985)  or  at  conferences  (cited 
above)  have  separately  described  only  gross  features  of  the  model  without  providing  necessary 
details.  Section  2  provides  the  basic  theory,  warped  coordinate  system,  mathematical 
formulation,  and  the  computational  algorithms  of  the  model.  The  second  objective  is  to  present 
the  results  of  a  critical  evaluation  of  the  model  using  the  MADONA  database  (section  3). 

Finally,  the  last  section  (section  4)  gives  conclusions  from  the  present  study  of  the  HRW 
model,  and  relates  these  results  to  a  companion  evaluation  of  the  model  in  a  report  by 
Williamson  et  al.  (2005). 


2.  Description  of  the  HRW  Model 


2,1  Gauss’s  Principle  of  Least  Constraint 


The  HRW  model  uses  Gauss’s  principle  of  least  constraint  directly  rather  than  the  more 
commonly  employed  Newton’s  laws  of  motion.  A  discussion  of  this  principle  applied  to  systems 
of  point  particles  can  be  found  in  Lanczos  (1966).  As  applied  to  non-viscous,  incompressible 
fluids,  the  principle  can  be  expressed  as 


8 


1  2 

\\\^p(a+s)  dz+\\pA  ' do 


=  0 


(1) 


where  £  denotes  a  variation  of  the  two  integrals  in  square  brackets.  The  vector  A  is  the 
acceleration,  dV  /  dt,  where  V  is  the  velocity  vector;  and  -g  is  the  acceleration  of  gravity,  the 
only  external  force  considered.  The  two  integrals  above  are  over  a  material  volume  (x)  and  its 
boundary  surface  (a),  respectively.  The  symbols  p  and  p  are  the  fluid  density  and  pressure 
respectively.  For  the  atmosphere,  p  and  p  are  related  by  the  equation  of  state  for  an  ideal  gas, 

P  =  pRT,  (2) 


where  R  is  the  specific  gas  constant  for  air,  and  T  is  the  absolute  temperature.  The  more 
conservative  potential  temperature  <9  is  related  to  T as 

e=T(pr,fipf”',  (3) 


where  pref  is  the  reference  pressure  taken  as  100  kPa  and  cp  is  the  specific  heat  at  constant 
pressure. 

Equating  the  variation  of  the  expression  in  brackets  to  zero  implies  that  the  air  motion  takes 
place  in  such  a  way  as  to  minimize  constraint  forces  arising  solely  from  kinematic  conditions.  In 
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this  case  the  constraints  are  the  boundary  conditions  and  conservation  of  mass,  which  for  an 
incompressible  fluid  is  expressed  as 

V-V=0.  (4) 

This  constraint  results  in  the  pressure  gradient  acting  as  a  force  of  constraint  necessary  to  enforce 
incompressibility.  Ball  and  Johnson  (1978)  discuss  this  result  in  greater  detail. 

In  order  to  account  for  atmospheric  thermal  effects,  additional  assumptions  are  needed.  First,  the 
variables  in  the  equation  of  state  are  decomposed  as 

P  =  Po+p\P  =  Po  +  P'J  =T0  +T'’&  =  00+0',  (5) 


where  the  variables  with  zero  subscripts  denote  ambient  or  mean  values  while  the  same  variables 
with  primes  represent  departures  from  the  ambient.  Secondly,  using  the  Boussinesq 
approximation,  equation  1  can  be  approximated  by 


S 


W\\p»(A+bg)  dT+\\p’A'dcj 


=  0 


(6) 


b  =  (p'/po)  =  (-d'/0o). 


(6a) 


where  -b  g  ,  the  effective  external  acceleration,  is  the  buoyancy  acceleration,  and  b  is  the 
buoyancy  parameter  defined  by  equation  6a.  A  detailed  discussion  of  the  Boussinesq 
approximation  is  provided  by,  e.g.,  Stull  (1988,  Chapter  3).  Likewise  the  effective  pressure  in 
equation  6  is  the  fluctuating  pressure  (//) ,  the  departure  from  the  ambient  pressure  (po).  It  is 
assumed  that  both  p  and  po  satisfy  the  hydrostatic  equation 

dp0  =  -p0g  dz  ,  dp'  =  -p’g  dz ,  (7) 


where  z  denotes  vertical  height.  Thirdly,  similar  to  equation  6  a  variational  expression  can  be 
applied  to  the  potential  temperature  field 


]_ 
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cW]2 

dt  y 


dz 


=  0, 


(8a) 


where  the  adiabatic  condition  has  been  assumed.  Equation  8a  implies  that 

o. 

dt  dt 


For  a  steady  state,  equation  8a  is  equivalent  to 


8 


=  0. 


(8b) 


(8c) 
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This  equation  is  needed  in  order  to  make  the  temperature  field  and  the  buoyancy  forces 
consistent  with  the  wind  field. 

For  the  application  of  the  basic  variational  principles,  equations  6  and  8c,  to  a  single  surface 
layer,  further  simplification  is  required.  This  single  surface  layer  is  of  constant  thickness  in  a 
direction  nonnal  to  the  terrain  surface  and  follows  the  warped  ground  surface.  The  second 
(surface)  integral  in  equation  6  contains  the  pressure  fluctuation  p' .  This  //  term  is  extremely 

difficult  to  estimate.  A  practical  alternative  is  to  reformulate  equation  6.  With  the  help  of 
Gauss’s  theorem,  equation  6  can  be  rewritten  for  the  single  surface  layer  as 

5  ilf{ \p*{As  +  An+bgs+bgn)  +V-jp'(4+4,)|jr  =0, 

where  the  subscripts  s  and  n  denote  the  surface  parallel  and  surface  normal  component, 
respectively.  The  integrand  in  the  above  equation  can  be  rewritten  as 

^a>  [As+^n+bss)  +(*>#„)  +2bA„-g„  +y-p'A„+y-p'As 

One  part  of  the  second-to-last  term  in  the  above  expression,  (V-p'An)  can  be  used  to  cancel  the 
last  tenn  within  the  brackets,  (p0bAn»gn)  since  dp'/ on  =  ~p'g„  as  implied  by  equation  7.  The 
absence  of  the  cross  product  of  A  n  and  bg  n  suggests  that  bg  n  may  be  ignored.  The  other  part 
of  the  second-to-last  term,  (p'V»An)  as  well  as  the  last  term,  (V-p'AJ  can  be  neglected  because 
p'  is  usually  small  Consequently  equation  6  is  further  simplified  as 

8  JIf(^  dx  =  0  (9) 

where  only  the  surface  parallel  component  of  buoyancy  bg  s  remains,  and  the  multiplicative 

constant  (po/2)  has  been  neglected.  The  surface  parallel  simplification  in  equation  9  assumes  that 
for  the  surface  layer  the  acceleration  normal  to  the  terrain  surface  is  not  affected  by  the  normal 
buoyancy  force  (bgn ) . 

2.2  Warped  Coordinate  System 

A  terrain-following,  warped,  non-orthogonal  coordinate  system  is  employed  in  the  HRW  model 
in  order  to  account  for  a  complex  and  varying  terrain  surface  as  accurately  as  possible.  Figure  1 
is  a  diagram  of  the  warped  coordinate  system.  The  upper  portion  of  figure  1  indicates  a 
Cartesian  coordinate  system  with  respect  to  a  horizontal  reference  plane  in  which  i ,  j  ,  and 
k  are  three  unit  vectors  along  the  three  orthogonal  directions  (x,  y,  and  z),  respectively.  Any 
position  on  the  terrain  surface  is  given  by  the  vector 

f  =  xi  +  y  j  +  h(x,  y)  k ,  (10) 


4 


where  h(x,  y)  is  the  surface  elevation. 


Figure  1.  The  terrain  following,  warped  and  non-orthogonal  coordinate 

system  with  three  base  vectors  (a, ,  a2 ,  and  a2 )  which  is  defined 

by  equations  11,  12,  and  13,  respectively.  V  is  a  wind  vector. 
A  computational  volume  element  with  top  (HMNP)  and  bottom 
(OEFG,  shaded)  is  also  shown,  see  text  for  details. 


The  lower  portion  of  figure  1  shows  three  base  vectors,  which  are  defined  by 

dr  -r  dh  |  _  |  r. 

al=  —  =  i+hxk,  hx=—,  ax=\ax\  =  ^l  +  hx 
ox  ox 


a2=^-  =  j  +  hyk  ,  hy=^~,  a2=\a2\  =  Juh 


dh 


dy 


dy 


a3  =  n  = 


h  =  ^j={-hxi -hj +  k\, 
V  a 


ci  —  1  +  h  ~  +  h  ,  I U-,  I  —  1  • 


'x  y  >  |w3 1 


(11) 

(12) 

(13) 


These  three  vectors  ax,a2,  and  d,  establish  a  warped  coordinate  system.  Notice  that  d]  and  a2 
are  parallel  to  the  surface,  but  are  not  unit  vectors.  The  vectors  d,  and  d2  are  not  perpendicular 
to  each  other,  and  hence,  the  warped  coordinate  system  is  not  an  orthogonal  coordinate  system. 
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The  vector  d,  is  a  unit  vector  and  is  perpendicular  to  the  terrain  surface  (d,  and  a2)  .  The  volume 
of  the  parallelepiped  formed  by  these  three  base  vectors  is 

a,  -[a2  x  a2  j  =  -Ja  .  (14) 

The  angle  between  d3  and  the  vertical  z  axis  is  a  as  indicated  in  figure  1, 

cosa  =  l/Va.  (15) 


It  is  useful  to  define  the  system  of  reciprocal  base  vectors  denoted  by  a  superscript 

d1  =  [a2  x  d3 )  /  \[a  =  — ^1  +  h2v  )/  -  hxhv  j  +  hx  k 


=  (d3  x  dj ) /  yfa  =  —  -hxhv  i  +  (l  +  h2x  +  hv  k 


d3  =  (dj  xd2)/  \la  =  /— 


-hx  i  -hvJ  +  k 


(16) 

(17) 

(18) 


The  reciprocal  vector  a 3  is  identical  with  d3  or  n  ,  the  unit  vector  nonnal  to  the  terrain  surface. 
The  reciprocal  vectors  d1  and  a 2  are  parallel  to  the  surface,  but  distinct  from  at  and  d2 .  Also 
d1  and  d2  are  not  perpendicular  to  each  other.  The  system  of  reciprocal  vectors  facilitates 
operations  because 

a,.  •  d '  =  Srs  r,s  =  1,2,3,  (19) 


where  8rs  =  1  if  s  =  r,  and  8rs  =  0  otherwise. 

A  three-dimensional  wind  velocity  vector  in  the  atmospheric  surface  layer  can  be  expressed  as 

V  =  V\  +  V %  +  V- %  =  V'a, .  (20) 


where  f'd,  =  OA  ,  V2a2  =  OC  ,  and  I/3d,  =  OD  as  illustrated  in  figure  1 .  The  vector  OB  is  the 

component  of  the  wind  velocity  along  the  terrain  surface.  The  sum  convention  for  repeated  index 
(i)  is  used  in  equation  20  and  hereafter.  Notice  that  (V)2  ^  (F1)2  +  (F2)2  +  (F3)2  due  to  non- 
orthogonality  of  the  coordinates.  The  component  F  is  called  the  surface  nonnal  velocity 
component  or  the  impaction  vertical  motion.  Note  that  F  is  not  the  vertical  velocity  (vv)  in 
terms  of  a  Cartesian  coordinate  system.  It  is  not  difficult  to  prove  that 


{(l  +  hl)u-hxhyv  +  hxw 

(21a) 

~hxhyu  +  (\  +  h2x)v  +  hyw 

(21b) 
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(21c) 


V3  =-i=[-hxu-hyv+w] 
v  a 

where  u  and  v  are  the  wind  velocity  components  along  x  and  y  direction,  respectively  (figure  1). 
The  surface  kinematic  boundary  condition  is  usually  expressed  by 

w(z  =  h)  =  V  ■  Vh(x,  y )  (22a) 


with  the  help  of  equation  21c,  the  above  condition  leads  to 

V3(z  =  h)  =  0  (22b) 

That  means  the  surface  kinematic  condition  imposes  zero  impaction  vertical  motion,  but  not  the 
Cartesian  vertical  velocity  (w)  at  the  terrain  surface.  As  shown  in  figure  1,  the  acceleration  of 
gravity  is 

~g  =  -gk  =  ~g  (cos  a  ■  n  +  sin  a  ■  s)  =  -  +  gs )  (23) 


where  -gn  and  —gs  denote  the  surface  normal  and  surface  parallel  component  of  -g , 

respectively.  The  term  S  is  a  unit  vector  along  the  OB  direction  in  figure  1 .  As  discussed 
earlier,  bgs  is  used  for  the  calculation  of  buoyancy  acceleration  in  equation  9.  Besides  the 

buoyancy  parameter  ( b ),  a  surface  heating  (dq)  is  introduced  in  the  model.  From  the  First  Law 
of  thennodynamics  for  an  ideal  gas,  the  heating  can  be  written  as 


dq  =  cp  dT  -  dp/p  =  cp  T(dQ/9) 


where  the  definition  of  0,  equation  3 ,  has  been  used,  see,  e.g.,  Wallace  and  Hobbs  (1977). 
Hence  for  the  surface  layer,  dq  can  be  related  to  the  buoyancy  parameter  ( b )  as  defined  by 
equation  6a 


dq  =  cPTsfc 


@sfc  6q 
^0  ' 


=  — c  T ,b . 

P  sfc  ' 


(24) 


where  the  subscript  “sfc”  denotes  the  surface  layer.  In  order  to  evaluate  the  integrals  of 
equations  8c  and  9,  we  use  a  computational  volume  element  for  the  single  surface  layer  as 
sketched  in  figure  1 .  This  volume  element,  also  called  local  flux  box,  consists  of  the  layer 
between  the  terrain  surface,  ( fi )  and  a  constant  normal  height  (d)  with  lateral  lengths  f  \  and  £ i 
as  indicated  in  figure  1 

f  j  =  OE  =  GF  =  la,  =  l -yjl  +  h;  (25a) 

£2=OG  =  EF  =  la2  =  l^\  +  h]  .  (25b) 
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where  /  =  OE'  =  E'F'  is  the  side  length  of  the  square  which  is  the  projection  of  the  volume 
element  on  the  horizontal  reference  plane.  The  volume  element  is  not  rectangular  and  its 
volume  (A  r)  is 


Ar  =  (/j  /aj)aj 


(/2  / a2)a2  xda2 


(26) 


For  reference,  the  lateral  faces  of  the  local  flux  box  are  numbered  counterclockwise  (1),  (2), 
(3),  and  (4)  as  shown  by  EFNM,  FGPN,  GOHP,  and  OEMH,  and  (5)  and  (6)  denote  top 
(HMNP)  and  bottom  (OEFG)  faces,  respectively  in  figure  1.  The  vector  expressions  of  the 
six  surface  areas  are 


aa 


AS\  =[(l2/ a2)a2xda2^=  Id  J~aal 
AS2  =  \jd -d2  x(/j /al)al =  Id  Jaa' 

AS 3  =  -[(4  I ai)^2  x^]3  =  -Id  \fcit 

AS4  =  -\jda2  x(/j /o1)a1]4  =  -Id  \faa~ 
AS 5  =  [(/j  /al)al  x(4  / a2 ) ^2 ]5  =  ^ 

AS  6  =  -[(/j  /  x(4  /  «2  )  «2  ]6  =  -l2  ^ 


J3 


J4 


(27a) 

(27b) 

(27c) 

(27d) 

(27e) 

(27f) 


where  a\a2,  and  a3  have  been  defined  by  equations  16,  17,  and  18  respectively.  The 
subscripts:  (1  through  6)  in  equation  27  mean  that  all  terrain  dependent  quantities  must  be 
evaluated  using  representative  (average)  values  on  that  surface.  The  surface  area  vector 
ASk  (k  =  1, . . .  6)  is  in  the  outward  normal  direction  of  each  face.  Using  the  divergence  theorem 
for  the  volume  element,  the  condition  of  mass  conservation  equation  4  can  be  approximated  by 


V -V 


-hW 


a ,  +  V^a2  +  V3a3 )  •  AA 


=  0. 


(28a) 


From  equations  19  and  27  the  above  condition  yields 


(28b) 


where  W is  only  a  notation  to  be  used  later.  The  term  (VaV3)6  in  equation  28b  vanishes  due  to 

equation  22b.  The  surface  normal  velocity  (V)  in  the  surface  layer  is  of  order  ( d/1 )  <0.1  with 

1  2 

respect  to  surface  parallel  components  ( V  and  V  term)  as  implied  by  equation  28.  Using 
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equations  19,  20,  and  27,  and  dyadic  (tensor)  analysis  (Tai,  1992,  for  example)  the  acceleration 
vector  for  the  local  flux  box  can  be  written  as 


A  =  v(ff)  =  V-(frJa,al)  =  -^'£\(f'V>3,aJ)AS' 

k  k 


With  the  help  of  equations  11,  12,  and  13  A  can  also  be  expressed  for  the  Cartesian  coordinate 
system  as 

2  =  [(vV'A\  +(r'v2-Ta\ -{r'v'4~a\ ~(r'v^\  +(r'w\\ 

+  ^]^V2Vlyfa^+(v2V2yfa)2-(v2Vlyfa^-(v2V2y[a^+(V2w)5]j 

+  ^=^{[V^(FIF\  +  F2F1/?J)]i  +[Va(F'F\+F2F2/jv)]2  (30) 

-  4a  (vlVxhx  +  V2Vlhv )]  -  [ sfa  (W2hx  +  V2V2hy )] 

+[w{v\  +  v\)]^, 

in  which  all  of  the  terms  of  V 3  V"1  ( m  =  1,  2,  or  3)  have  been  neglected  since  V3  «  V1  or  V2  and 
equation  28b  has  been  used.  V  can  be  large  only  if  hx  and/or  hv  in  equation  2 lc  is  large. 

2.3  Empirical  Surface  Layer  Profiles 

Accurate  calculation  of  the  integrals  in  equations  8c  and  9  for  a  volume  element  requires  use  of 
surface  nonnal  profiles  of  temperature  and  wind  velocity  in  the  surface  layer.  The  HRW  model 
assumes  that  the  following  empirical  profiles  in  the  z  direction  may  still  be  applied  locally  in  the 
n  direction  (figure  1).  The  model  has  adopted  profiles  for  wind  speed  and  potential  temperature 
as  follows: 
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V(z)  =  Vr 


f  z^ 


\ZrJ 


(31) 


and 


6(z)  =  6r  +  y(z  -  zr)  , 


(32) 


where  Vr  and  6,  denote  the  wind  speed  and  the  potential  temperature  at  a  reference  height  (zr) 
which  can  be  the  top  of  the  volume  element.  From  these  equations  note  that  the  potential 
temperature  is  assumed  linear  with  height,  /being  the  gradient,  and  that  the  wind  speed  follows  a 
power  law  (e.g.,  Arya,  1999).  The  value  of  the  exponent  m  can  be  estimated  from  surface  layer 
similarity  theory,  which  defines  a  non-dimensional  wind  shear 


f  Z  ^ 
\Lj 


kz  dV 
ut  dz 


(33) 


Hence 


V(z)  = 


LA 

L) 


Vn 


f  Z  ^ 
\Lj 


(34) 


where  z0  is  the  roughness  length,  L  is  the  Obukhov  length,  and  ij/m  is  a  newly  defined  function; 
see  Garratt  (1992)  for  details.  Panofsky  and  Dutton  (1984)  have  shown 


m  = 


6 

t  m 


f  Z  ^ 


ln  — 

zn 


f  z  ^ 

\Lj 


(35) 


The  empirical  wind  shear  and  wind  speed  profiles  of  Businger  et  al.  (1971)  have  been  used  by 
the  HRW  model.  It  is  impractical,  however,  to  use  their  original  formulation  since  L  and  hence 


the  non-dimensional  height  £ 


f  z  ^ 

"  L 


are  usually  unknown.  Instead,  an  empirical  expression  for 


<j)m  in  tenns  of  the  Richardson  number  (Ri)  can  be  used.  For  example,  Ball  and  Johnson  (1978) 
used  the  following  expression: 


^  =(1-12  Rip;  Ri  <  0.0 

=  (l-37?i)  1 ;  0.0<R/<0. 03571 
=  0.88(1-67?/)  1 ;  0.03571  <  Ri  <  0.1246 
=  1.75(1-47?/)  1 ;  0.1246  <7?/  <  0.25 


(36) 
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in  which 


Ri  =  (Ri)b  +  (Ri)c, 


where  the  bulk  Richardson  number,  (Ri)b,  is 


mb 


'g)<Oo-8) 

U  J  n2 


(37) 


(38) 


(Ri)c  in  equation  37  represents  the  Richardson  number  for  curved  flows.  This  additional 
Richardson  number  was  originally  suggested  by  Bradshaw  (1969)  for  curvature  effects 
of  curved  flows, 


(Ri)c=2S(l  +  S), 


(  dv\ 

UJ 

KdR  j 

(39) 


Here  R  is  the  radius  of  streamline  curvature.  If  the  average  surface  nonnal  acceleration  A„  in  the 
layer  is  used  and  is  assumed  to  be  purely  centripetal  (-VV/R),  and  the  average  value  of  V(dVldz) 
of  equation  3 1  is  used,  then  S  can  be  approximated  by 


5  = 


2  Az. 


V: 


2.4  Computational  Aspects  and  Algorithm 

The  computational  aspects  of  the  HRW  model  are  illustrated  in  figure  2,  a  simplified 
computational  flow  chart.  As  indicated  in  figure  2  there  are  three  main  stages:  (1)  data  input 
and  model  setup;  (2)  relaxation  scheme  through  direct  variation;  and  (3)  model  data  processing, 
output  and  plotting  programs.  These  primary  stages  are  discussed  in  the  following  sections. 
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Stage  1 

Data  Input  & 
Model  Setup 


Algorithm 

Major  routines 

1.  Input  data 

Read_Hinfil,  Read_Terrain, 
Read_Morpholcgy 

2.  Model  setup 

Inputs,  Setup 

T 


Algorithm 

Major  routines 

1.  Obtail  local  quantities 

Locqua 

2.  Compute  integrals  in 
(9)  and  (8c) 

Accsqr,  Tpflux 

3.  Apply  steepest  descent 

corrections  to  wind  and 
temperature  field 

Delta 

Stage  2 
Relaxation 
Scheme 


Stage  3 

Processing  & 
Data  Output 


Algorithm 

Major  routines 

1.  Calculate 
turbulence  quantities 

Turb 

2.  Output  model  data 

Output,  Trbout 

Figure  2.  Computational  flow  chart  for  the  HRW  model.  This  chart 

illustrates  the  three  essential  stages  and  their  major  routines. 


2.4.1  Data  Input  and  Model  Setup 

To  initiate  a  model  run,  necessary  input  data  are  required  to  provide  control  and  driving 
parameters.  Table  1  (upper  part)  lists  the  required  input  data.  First,  there  are  a  number  of 
control  parameters  which  include  the  setting  of  the  computation  domain,  grid  system  and 
resolutions,  important  constants,  data  handling  and  storage,  etc.  Secondly,  the  computation 
requires  input  of  digitized  terrain  data  at  every  grid  point,  h(i,j)  i,j  =  1,  2,  3,  -  - ,  N.  The  input 
file  of  terrain  height  for  the  major  routine  (Read-Terrain,  in  figure  2)  should  be  properly 
formatted.  The  additional  vegetation  data,  i.e.,  the  vegetation  element  heights,  hveg(i,j ),  are  also 
needed.  In  such  a  case,  the  effective  terrain  heights,  heff(i,j),  will  be  calculated  as 

Kff  (b  j)  =  KU  j)  +  0.70  *  hveg O',  j). 


An  empirical  expression  for  the  roughness  length  (zq)  over  vegetation  has  been  used 


z0(i,  j)  =  0.014  +  0.14*  hveg  0,  j) , 


where  z0  is  in  meters,  see  Cionco  (1983). 
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Table  1 .  A  list  of  input  and  output  data  for  the  HRW  model. 


Input 

Category 

Notation  (unit) 

Remark 

1 

Control  parameters 

D(x,  y),  Computation  domain 

And  flags 

Ax,  Ay  resolutions 
grid  system,  relaxation 
processes,  etc 

typically  Ax  =  A y=  100  m 

2 

Digitized  terrain  data 

Khj),  (m) 

II 

0 

K> 

II 

0 

3 

Vegetation  data 

hveg  0*,  /),  (m) 

Vegetation  element  height 

4 

Surface  station  observation 

U  (at  10  m),  (ms1) 

Constant  value  within  D 

Tsfc  (at  10  m),  ( K) 

P( 3),  7X3),  h( 3) 

pressure,  temperature  and 

5 

Upper  air  sounding 

P( 2),  7X2),  h( 2) 

geopotential  height  at 

Ml),  7X1),  h( l) 

three  levels,  e.g.,  at  700,  850,  mb 
and  at  the  surface 

Output 

1 

Calculated  wind  field 

v  a,  j)  (ms  ') 

defined  by  (21) 

2 

Calculated  potential 
temperature  field 

dsfcd,  j)  (K) 

calculated  by  (24b) 

3 

Richardson  number 

RKU  j) 

defined  by  (37) 

4 

Friction  velocity 

u,(Uj)(ms~') 

calculated  by  (34) 

5 

Wind  profile  exponent 

m 

calculated  by  (35) 

As  indicated  in  table  1 ,  a  minimal  set  of  meteorological  data  is  also  required  to  initialize  the 
model.  The  initial  wind  field  (Vsfc)  is  obtained  either  by  local  meteorological  observation  of 

wind  speed  and  direction  (at  10  m  height)  or  by  an  estimation  of  the  geostrophic  wind  near  the 
surface.  The  model  assumes  a  homogeneous  initial  wind  field,  i.e.,  the  initial  velocity  V, ,  is 

constant  throughout  the  computational  domain.  On  the  other  hand,  the  initial  potential 
temperature,  Osfc(i,j),  varies  within  the  domain.  To  calculate  the  initial  field  of  6(i,j),  the 
following  algorithm  has  been  adopted. 

(a)  Using  the  upper  air  sounding  as  the  input  data,  the  ambient  potential  temperature  ( do)  as 
defined  by  equation  5  is  derived  by  extrapolating  the  background  potential  temperature  profile  to 
the  surface.  For  example 

K  =  a™  +{h* (4°) 

"850  "700 

where  6^50  and  6(00  are  the  potential  temperature  at  the  850  and  700  mb,  respectively.  Further, 
hsfc,  /?X50,  and  him  refer  to  the  geopotential  height  at  the  surface,  850  and  700  mb,  respectively, 
and  hsfc  is  also  the  altitude  of  the  sounding  station. 
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(b)  Also  using  the  upper  air  sounding  and  the  terrain  height  data,  h(i,j),  the  surface  pressure  at 
each  grid  point  is  calculated  by  a  linear  extrapolation,  e.g., 


P(U  j)  =  850  +  [h(i,j)  -  h&50  ]  ^5° 


h  -h 

"850  "700 


where  the  p(i,j)  has  the  unit  of  mb. 


(c)  To  calculate  the  surface  heating  (A Q)  as  an  input  parameter  from  equation  24 


(24a) 


where  the  surface  temperature  ( Tsfc)  and  the  surface  potential  temperature  ( 6sfc)  are  obtained  from 
observation  of  a  surface  meteorological  station. 

(d)  It  is  assumed  that  the  value  of  A Q  from  a  single  (station)  observation  is  representative  for  the 
whole  domain.  This  assumption  means  the  right  hand  side  of  equation  24a  can  be  applied  to  any 
grid  point  (/,/).  Thus 


^(bj)  =  #0(bi)  +  A£> 


0o(Uj) 


0()(i,j)  +  AQ 


10.286 


1000 


p(i,j) 


(24b) 


where  the  ambient  background  potential  temperature  at  a  gridpoint  (/,})  can  be  extrapolated  from 
the  sounding  data,  equation  40. 


ojij) = 4, +p(-.y)-*,„]  f”  _y 

"850  "700 


(24c) 


Lanicci  (1985)  has  discussed  the  calculation  of  the  A  Q  and  has  suggested  some  modifications. 
For  a  location  where  hsfc  >  A 8 so,  we  will  use  the  sounding  data  at  higher  levels,  e.g.,  500  and 
700  mb  instead  of  at  lower  levels,  e.g.,  700  and  850  mb. 

2.4.2  Relaxation  Scheme 

As  indicated  by  figure  2,  the  central  stage  of  the  HRW  model  computation  is  the  relaxation 
scheme,  which  is  the  workhorse  of  the  model.  The  scheme  implements  the  calculation  of  the 
integrals,  equations  8c  and  9  as  well  as  the  variational  technique  to  find  the  minimum  of 
equation  9.  As  mentioned  in  section  2.2,  the  evaluation  of  the  integral  in  equation  9  adopts 
computational  elements  called  flux  boxes.  Consequently,  the  total  constraint  integral  of 
equation  9  can  be  expressed  as  a  sum  over  all  flux  boxes  in  the  modeled  area 

*r=2X=ZC+(«.),2  •  <9a> 

hj  Uj  l’J 
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For  variational  purposes  we  consider  the  Rfj,  called  residues,  as  functions  of  only  the  two 
velocity  components  ( V  ,  V  )ij  at  the  grid  point  (/,/).  Denoting  the  two  velocity  components  at 
each  grid  point  by  V^ik  =  1,2),  the  unit  vector  of  global  steepest  descent  in  this  two- 

dimensional  velocity  space  is  written  as 


z 


dKj 


\_  • 
2 


(41) 


The  steepest  descent  method  computes  a  velocity  correction  (Ah). )  proportional  to  the  steepest 

descent  unit  vector  for  each  relaxation  sweep  over  the  entire  grid.  This  velocity  correction  has 
been  chosen  as 

AFj  =  Cr'KylN-Mnl 


where  Cr  is  a  selected  fraction  of  a  velocity  scale,  Vr.  Currently  C,  =0.02  and  Vr  is  the  constant 
initial  wind  speed  or  2  ms-1,  whichever  is  greater.  N and  Mare  the  grid  dimensions,  currently 
N  =  M  =  5  \ .  All  velocity  corrections  are  applied  simultaneously  at  the  end  of  each  sweep. 
This  implies  the  velocity  field  after  each  relaxation  step  is 


Kj  =  Kj+aKp 


Vij  =  +AVtj. 


This  relaxation  algorithm  apportions  larger  corrections  to  regions  of  the  grid  where  the  constraint 
integral  is  most  sensitive  to  changes  in  the  velocity  field.  The  relaxation  scheme  calculates  the 
total  constraint  Rj  at  each  relaxation  sweep  and  saves  the  wind  field  for  the  minimum  value  of 
Rt  achieved.  Each  relaxation  sweep  also  adjusts  the  potential  temperature  field  by  use  of 
equation  8c.  Consequently  each  relaxation  sweep  yields  altered  potential  temperature  and 
velocity  fields  until  equation  9  is  satisfied. 

2.4.3  Data  Output 

The  final  stage  of  the  model  computation  produces  various  computational  results.  As  shown  in 
table  1,  the  output  data  comprises  at  least  five  categories.  The  most  important  output  is,  of 
course,  the  surface  layer  wind  field.  It  should  be  stressed  that  the  calculated  windfield  refers  to 
the  three  components  in  the  terrain  following  warped  coordinate  system  since  the  model  is 
formulated  according  to  such  a  unique  coordinate  system.  It  is,  however,  easy  to  transform  these 
three  components  ( V  ,  V  ,  V  )  into  Cartesian  wind  components  (u,  v,  w)  as  indicated  by 
equation  2 1 .  The  second  category  of  output  is  the  calculated  field  of  surface  potential 
temperature  which  can  also  be  transformed  into  surface  temperature  easily.  Finally  the  other 
three  categories  of  output  are  related  to  turbulence  characteristics  ( ut  and  Ri)  and  the  exponent 
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(m)  of  the  power  law  of  wind  profiles  in  the  surface  layer.  During  the  past  decade  several 
versions  of  algorithms  and  programs  for  output  displays  and  plotting  have  been  tested  and 
established.  Those  aspects  related  to  the  software  developments  will  not  be  discussed 
in  this  report. 


3.  Results  of  the  Model  Evaluation 


3.1  The  MADONA  Data 

The  micrometeorological  data  used  in  this  evaluation  study  are  based  on  measurements  during 
the  field  study  of  MADONA  in  September  1992.  Cionco  et  al.  (1999)  have  provided  a 
comprehensive  overview  of  the  MADONA  project  including  its  database.  The  MADONA  field 
site  is  located  at  Porton  Down,  near  Salisbury,  United  Kingdom.  From  14  to  23  September  1992, 
there  were  ten  days  of  intensive  measurements.  Extensive  high-resolution  micrometeorological 
data  were  collected  to  complement  the  diffusion  data.  Fifteen  sets  of  standard  wind  speed  (three 
cup  anemometer)  and  wind  direction  (vane  anemometer)  sensors  were  deployed  at  14  locations 
by  the  UK  Chemical  and  Biological  Defence  Establishment  (CBDE,  now  part  of  the  Defence 
Evaluation  and  Research  Agency).  These  wind  sensors  were  calibrated  before  and  after  the  field 
measurements.  Figure  3  illustrates  the  locations  of  these  wind  measurement  stations  (Ml  through 
Ml 5)  plotted  on  terrain  contours.  All  wind  sensors  except  M9  were  mounted  at  10  m  above  the 
ground.  The  M9  station  was  installed  at  30  m  above  the  ground  at  the  same  tower  as  M8  to 
provide  dual  measurements  at  the  same  location.  The  area  covered  by  figure  3  is  5  by  5  km, 
which  is  the  domain  for  the  HRW  simulations.  The  wind  data  from  Ml  and  M3  are  excluded  for 
model  evaluation  because  their  locations  are  very  close  to  the  model  boundary,  and  the  data  from 
M8/9  are  not  used  due  to  technical  reasons.  Altogether  the  measured  wind  data  from  eleven 
locations  (Station  M2,  M4-7,  and  Ml 0-1 5)  are  used  for  the  comparison  between  the 
measurements  and  the  HRW  simulations.  All  the  wind  data,  as  well  as  other  MADONA  data 
including  the  area  terrain  data,  have  been  officially  documented  on  CD-ROM.  The  CD-ROM 
includes  the  averages  of  wind  speed  and  direction  for  5  min,  10  min,  30  min,  and  60  min,  among 
other  statistics.  The  5  min  averages  of  wind  measurements  at  10  m  above  the  ground  from  these 
1 1  locations  are  those  used  for  the  present  study.  Unfortunately,  there  were  no  corresponding 
temperature  sensors  at  the  same  1 1  stations  provided  by  CBDE  due  to  financial  constraints. 

The  CD-ROM  is  available  on  request  from  the  RISO  National  Laboratory,  P.O.  Box  49, 

Roskilde,  Denmark. 


16 


Figure  3.  Locations  of  MADONA  measurements  sites  (Ml  through  M15)  plotted 
on  terrain  contours. 

As  shown  in  figure  3,  the  topography  at  the  MADONA  site  is  gently  rolling  terrain  with  a  ridge 
running  southwest  to  northeast  (approximately  230°-50°)  with  higher  terrain  at  each  end.  Total 
relief  for  the  domain  is  about  100  m  from  80  m  northwest  to  180  m  near  Tower  Hill.  In  such  a 
small  area,  however,  the  wind  distribution  inhomogeneity  is  significant,  as  shown  by  Cionco  et 
al.  (1999,  their  figures  7  and  8).  Generally,  the  wind  speeds  over  the  ridge  (e.g.,  M14)  are 
stronger,  while  over  the  lower  land  (e.g.,  Ml 3)  they  are  weaker. 

3.2  The  Model  Simulations 

During  the  ten  days  of  intensive  measurements  from  14  to  23  September  1992,  a  series  of 
boundary  layer  soundings  were  launched  for  each  date  from  the  middle  of  the  MADONA 
experimental  area,  as  described  by  Cionco  et  al.  (1999).  Also,  traditional  upper  air  soundings 
from  the  UK  Meteorological  Office  Larkhill  station  and  surface  meteorological  observations 
from  CBDE’s  Met  office  were  available  during  the  experiment.  These  sounding  and  surface 
observation  data  (surface  temperature  at  10  m  AGL  in  particular)  provide  the  necessary  inputs 
for  the  HRW  simulations.  The  wind  data  measured  at  Station  M10  are  used  to  initialize  the  wind 
field  for  model  runs.  Table  2  lists  the  39  model  runs  chronologically  during  these  10  days.  The 
surface  heating  (AQ)  in  table  2  is  calculated  from  equation  24a.  The  wind  speed  and  direction  in 
table  2  are  the  initialized  values  which  are  the  5  min  averages  of  the  observed  wind  at  Station 
M10,  located  in  the  middle  of  the  diffusion  study  area. 


17 


Table  2.  The  39  cases  for  the  HRW  simulations.  A Q  is  defined  by  equation  24a  in  the  text. 


CaseNo. 

Day 

(1992) 

Local 

Time 

Afi 

(K) 

Wind  Speed 
(ms-1) 

WindDirection 

(degree) 

1 

14  Sept 

14:10' 

9.21 

5.4 

261 

2 

15:55' 

9.78 

7.2 

280 

3 

17:20’ 

2.58 

4.0 

270 

4 

15  Sept 

12:10' 

-5.10 

8.0 

240 

5 

13:55’ 

-2.50 

7.5 

240 

6 

15:55' 

-2.70 

7.4 

245 

7 

17:40’ 

2.98 

5.5 

240 

8 

16  Sept 

10:30’ 

-2.45 

3.0 

120 

9 

13:55’ 

-3.45 

0.3 

90 

10 

16:10’ 

-5.45 

3.0 

160 

11 

17:50' 

-5.52 

2.0 

160 

12 

17  Sept 

10:30’ 

-6.54 

5.0 

130 

13 

14:15' 

-5.65 

1.0 

180 

14 

17:25' 

-4.02 

5.5 

180 

15 

18  Sept. 

10:30' 

-0.54 

4.5 

280 

16 

13:10' 

-0.88 

4.0 

255 

17 

14:55' 

-2.34 

3.0 

250 

18 

18:15' 

0.00 

2.5 

170 

19 

19  Sept 

10:30’ 

3.88 

4.0 

260 

20 

12:10' 

6.12 

5.0 

270 

21 

14:00' 

6.59 

5.0 

270 

22 

16:10’ 

5.94 

4.5 

275 

23 

18:25' 

3.56 

1.5 

250 

24 

20  Sept 

10:30’ 

3.73 

4.5 

220 

25 

12:40' 

10.26 

2.5 

225 

26 

14:45' 

5.39 

2.5 

225 

27 

15:55' 

3.86 

2.0 

200 

28 

21  Sept 

10.30' 

-0.20 

7.5 

70 

29 

14:25' 

1.29 

4.5 

110 

30 

15:25' 

-1.35 

3.5 

140 

31 

18:00’ 

-1.83 

3.5 

160 

32 

22  Sept 

10:30' 

1.11 

0.5 

320 

33 

13:30’ 

7.32 

1.2 

295 

34 

15:10' 

2.20 

2.5 

360 

35 

23  Sept 

10:30' 

5.44 

6.0 

280 

36 

12:00’ 

8.78 

7.5 

275 

37 

14:15' 

8.76 

8.5 

270 

38 

16:10’ 

9.71 

5.0 

280 

39 

17:40’ 

1.48 

2.5 

250 
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As  mentioned  before,  the  domain  of  the  model  simulation  is  5  by  5  km  with  a  spatial  resolution 
of  100  m.  Thus,  there  are  2601  (5 1  times  51)  grid  points  covered  by  the  computation  domain. 
Only  the  simulated  wind  speeds  and  directions  at  the  1 1  gridpoints  corresponding  to  the  1 1 
measurement  locations  have  been  chosen  for  comparison.  To  evaluate  the  model  performance, 
the  simple  linear  regression  analysis  is  used  as  follows 

y  =  ax  +  b  (42) 

where  x  is  the  observed  value  of  either  wind  speed  or  wind  direction,  and  y  is  the  corresponding 
simulated  value  from  the  model  runs.  The  terms  a  and  b  represent  the  slope  and  intercept  of  the 
regression  line,  equation  42,  respectively.  For  each  run  (case)  listed  in  table  2,  there  are 
generally  1 1  pairs  (datapoints)  of  (x,  v).  However,  there  are  only  10  datapoints  for  case  9,  and  9 
datapoints  for  case  23.  In  these  two  cases,  the  observed  wind  speeds  from  M12  (case  9)  and 
from  both  M12  and  M13  (case  23)  were  zero.  Consequently,  the  wind  direction  is  undetermined 
for  these  two  cases.  Therefore,  the  total  number  of  (x,  y)  datapoints  from  the  39  cases  is  426  for 
the  simultaneous  wind  direction  and  speed  evaluation.  For  a  single  station,  except  M 12  and 
Ml 3,  on  the  other  hand,  there  are  39  datapoints  for  the  linear  regression  analysis.  A  correlation 
coefficient  (r)  and  a  standard  deviation  (S)  for  a  set  of  (x,  y )  datapoints  can  also  be  calculated. 

The  linear  regression  should  yield  a  =  1,  b  =  0,  r  =  1,  and  S  =  0  for  an  ideal  or  perfect  model.  In 
reality,  however,  no  model  is  perfect.  The  values  of  a,  b,  r,  and  S  can  help  evaluate  the  model 
performance  under  various  circumstances. 

3.3  The  Results  of  Linear  Regression  Analysis 

(1)  Figures  4a  and  4b  present  the  simulated  winds  from  the  39  cases  versus  the  observed  winds 
in  scatter  diagrams  with  the  regression  line  equation  42  superimposed.  Figure  4a  is  for  the  wind 
direction  comparison,  while  figure  4b  is  for  the  wind  speed  comparison.  Both  figures  show  the 
plot  for  all  1 1  stations  together  as  well  as  for  individual  stations.  Table  3  lists  the  values  of  a,  b, 
r,  and  S  corresponding  to  figures  4a  and  4b.  From  the  total  426  datapoints,  the  values  of  a,  b,  r, 
and  S  are  1.004,  2.000,  0.970,  and  16.3°,  respectively  for  the  wind  direction,  and  0.868,  0.902, 
0.857,  and  1.02  ms-1,  respectively  for  the  wind  speed,  as  shown  in  table  3.  This  result  has  been 
partially  reported  by  Cionco  and  Byers  (1995).  Therefore,  the  HRW  model  simulations  as  a 
whole  appear  to  provide  quite  high  positive  correlations,  especially  for  the  wind  direction.  For 
individual  stations,  r  is  always  greater  than  0.96  for  wind  direction  and  greater  than  0.80  for 
wind  speed,  as  indicated  in  table  3.  The  model  performs  fairly  well  for  each  of  the  1 1  locations, 
although  it  appears  to  perfonn  better  for  certain  locations  (e.g.,  M10,  Ml  1)  than  for  other 
locations  (e.g.,  M15,  M14),  as  seen  from  figure  4  and  table  3.  The  model  does  not  demonstrate 
significant  failure  at  any  particular  location.  It  is  also  evident  that  the  model  simulations  appear 
better  for  wind  direction  than  for  wind  speed.  As  for  an  individual  case,  the  model  sometimes 
simulates  a  case  satisfactorily,  but  not  always.  A  statistical  result  for  each  individual  case  is 
meaningless  due  to  small  number  of  data  points,  and  hence  is  not  presented. 
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Figure  4a. 


Scatter  diagrams  for  simulated  wind  direction  (y,  degree)  versus  observed  values  (x,  degree). 
The  straight  lines  result  from  the  linear  regression  analysis,  equation  42. 
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Figure  4b.  Scatter  diagrams  for  simulated  wind  speed  (y,  ms  '  )  versus  observed  values  (x,  ms  '). 
The  straight  lines  result  from  the  linear  regression  analysis,  equation  42. 
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Table  3.  The  value  of  n  (data  pair  number),  a  (slope),  b  (intercept),  r  (correlation  coefficient),  and  S  (standard 
deviation)  from  the  linear  regression  equation  42,  for  wind  direction  and  wind  speed,  respectively. 


Station 

n 

Wind  direction  (degree) 

Wind  speed  (ms  ') 

a 

b 

r 

5” 

a 

b 

r 

5” 

M2 

39 

1.0116 

-5.3141 

0.9716 

15.9 

0.8169 

0.9626 

0.8739 

1.02 

M4 

39 

1.0288 

-7.0611 

0.9762 

15.1 

0.7601 

1.0301 

0.8358 

1.00 

M5 

39 

1.0241 

-8.4746 

0.9759 

14.7 

0.7964 

0.8659 

0.9397 

0.57 

M6 

39 

1.0795 

-11.2393 

0.9676 

17.9 

0.8728 

0.8877 

0.8582 

1.07 

M7 

39 

1.0278 

5.9955 

0.9697 

17.0 

0.8448 

1.2949 

0.8942 

0.92 

M10 

39 

1.0386 

-10.035 

0.9926 

8.3 

0.8951 

0.5111 

0.9586 

0.52 

Mil 

39 

0.9833 

7.5243 

0.9862 

11.0 

0.3993 

0.9655 

0.9255 

0.73 

M12 

37 

1.0344 

-5.6589 

0.9778 

14.8 

0.8725 

0.6401 

0.8629 

0.90 

M13 

38 

0.9391 

20.405 

0.9742 

14.9 

1.0292 

1.0245 

0.8502 

1.07 

M14 

39 

0.9631 

15.7441 

0.9775 

14.5 

0.8834 

0.8918 

0.8038 

1.40 

M15 

39 

0.9255 

7.0691 

0.9613 

18.7 

1.2376 

-0.5478 

0.8708 

1.11 

All 

426 

1.0043 

2.0001 

0.9701 

16.3 

0.8681 

0.9024 

0.8567 

1.02 

(2)  To  further  investigate  the  model  perfonnance,  comparisons  of  wind  field  between 
simulations  and  observations  for  different  wind  directions  and  different  wind  speeds  have  been 
conducted.  The  39  cases  listed  in  table  2  have  been  grouped  into  three  categories  of  wind 
direction  and  three  categories  of  wind  speed.  The  criteria  for  wind  direction  categorization  take 
into  account  the  general  ridge  direction  in  the  MADONA  area,  which  is  approximately  230°- 
50°,  as  mentioned  earlier  (see  figure  3).  The  three  wind  direction  (WD)  categories  are  defined  as 

1)  Parallel  to  the  ridge  (14  cases), 

200°  <WD<  260°  or  20°  <  WD  <  80°  (43a) 

2)  Perpendicular  to  the  ridge  (10  cases) 

1 10°  <  WD  <  170°  or  290°  <  WD  <  350°  (43b) 

3)  Slant  to  the  ridge  (15  cases) 

350°  <  WD  <  20°  or  80°  <  WD  <  1 10° 
or  170°  <  WD  <  200°  or  260°  <  WD  <  290°  (43c) 

The  three  wind  speed  ( WS)  categories  are  defined  more  or  less  subjectively.  They  are 

1)  Light  wind  (15  cases), 

WS  <  3.0  ms-1  (43d) 

2)  Moderate  wind  (16  cases), 

3.0  ms-1  <  WS  <  6.0  ms-1  (43e) 
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3)  Strong  wind  (8  cases), 


WS  >  6.0  ms-1  (43f) 

The  results  of  these  comparisons  are  presented  in  table  4.  They  are  also  shown  in  figure  5  for  the 
three  wind  direction  categories,  and  in  figure  6  for  the  three  wind  speed  categories,  respectively. 


Table  4.  Same  as  table  3  except  for  three  wind  direction  categories  as  defined  by  equation  43  (a,  b,  c)  and  three 
wind  speed  categories  as  defined  by  equation  43  (d,  e,  f). 


Category 

n 

Wind  Direction 
(degrees) 

Wind  Speed 
(ms-1) 

a 

b 

r 

5” 

a 

b 

r 

5” 

Wind  dir  ection: 
Parallel 

152 

0.9684 

11.7935 

0.9335 

17.8 

0.7946 

0.9584 

0.8878 

0.90 

Perpendicular 

110 

1.0028 

-1.5884 

0.9737 

15.4 

0.7447 

1.1904 

0.7117 

0.83 

Slant 

164 

0.9530 

16.7129 

0.9484 

14.9 

0.9111 

1.0255 

0.8144 

1.12 

Wind  speed: 

Light 

162 

0.9787 

6.5416 

0.9452 

21.8 

0.7949 

0.9477 

0.6617 

0.79 

Moderate 

176 

1.0271 

-1.5462 

0.9829 

12.4 

0.8645 

1.0031 

0.7293 

1.04 

Strong 

88 

1.0013 

1.1356 

0.9888 

10.0 

0.1479 

6.3077 

0.1507 

1.04 
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Figure  5.  Same  as  figure  4a  except  for  three  wind  direction  categories  as  defined  by  equation  43a,  b,  c. 
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Figure  6.  Same  as  figure  4b  except  for  three-wind  speed  categories  as  defined,  by  equation  43d,  e,  f. 
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It  appears  from  table  4  and  figure  5  that  the  model  perfonned  better  for  the  perpendicular  wind 
direction  category  with  the  highest  value  of  r  (0.974)  for  wind  direction  and  the  lowest  value  of  S 
(0.83  ms”1)  for  wind  speed.  It  is  conceivable  that  the  airflow  will  experience  more  acceleration  or 
deceleration  when  its  direction  is  perpendicular  to  a  ridge.  Consequently,  the  magnitude  related 
to  the  A  term  in  equation  9  will  be  more  significant  as  compared  to  the  parallel  or  slant  wind 

directions.  For  the  three  categories  based  on  wind  speed,  figure  6  and  table  4  reveal  an 
unsatisfactory  situation  for  the  strong  wind  category.  The  scatter  diagram  at  the  right  lower 
comer  of  figure  6  shows  that  the  datapoints  are  scattered  around  the  regression  line  with  low 
value  of  r  (0.1507)  and  high  value  of  S  ( 1 .04  ms”1)  although  the  simulations  for  the  wind 
direction  appear  much  better  than  for  the  wind  speed.  This  result  receives  more  discussion  in  the 
following  subsection. 

(3)  Finally,  the  39  cases  have  also  been  grouped  in  terms  of  atmospheric  stability  for  model 
performance  evaluation.  Based  on  the  value  of  A Q,  equation  24a,  the  following  three  stability 
categories  are  defined  as 


1)  Unstable  (19  cases), 


AQ  >  2.0° 


2)  Near-neutral  (9  cases), 

|A0|<  2.0° 


(44a) 

(44b) 


3)  Stable  (11  cases), 


AQ  <  -2.0° 


(44c) 


The  results  for  the  three  stability  categories  are  shown'in  figure  7  and  table  5.  Among  the  three 
atmospheric  stability  categories,  the  model  appears  to  perfonn  the  best  for  the  near-neutral 
stratification  as  compared  to  unstable  or  stable  stratification.  For  example,  both  values  of  r  for 
wind  direction  and  speed  are  the  highest  for  the  near-neutral  category.  This  result  will  be 
discussed  in  the  next  section. 
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Figure  7.  Same  as  figure  4  except  for  three  stability  categories  as  defined  by  equation  44a,  b,  c. 
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Table  5.  Same  as  table  3  except  for  three  atmospheric  stability  categories  defined  by  equation  44a,  b,  c. 


Category 

11 

Wind  direction  (degree) 

Wind  speed  (ms  ') 

a 

b 

r 

5” 

a 

b 

r 

5” 

Unstable 

207 

0.8949 

33.8784 

0.8881 

15.7 

0.9291 

0.9063 

0.8561 

1.00 

Near-neutral 

99 

1.0112 

-0.4089 

0.9759 

16.4 

0.7165 

0.8671 

0.8711 

0.67 

Stable 

120 

0.9379 

11.1184 

0.9444 

16.0 

0.7445 

1.4706 

0.8512 

1.01 

3.4  Discussion 


(1)  As  presented  in  the  last  subsection,  the  HRW  model  simulates  the  wind  field  better  for  the 
near-neutral  condition  than  for  either  stable  or  unstable  conditions.  This  is  because  under  neutral 
conditions  ( b  =  0)  the  buoyancy  acceleration  plays  no  role  and  the  wind  field  is  influenced  by 
topography  only.  This  result  implies  that  the  buoyancy  (b)  term  is  important  and  that  a  realistic 
representation  of  this  term  for  the  model  is  a  challenge.  As  defined  by  equations  24  and  24a,  the 
calculation  of  the  buoyancy  factor  (b)  or  the  surface  heating  (A0  involves  the  background 
(reference)  potential  temperature  (0o),  the  surface  air  temperature  (Tsfc),  and  the  related  surface 
air  potential  temperature  {6sfc).  The  calculation  of  these  three  input  parameters  is  not  trivial. 
Currently  0O  is  calculated  through  a  linear  extrapolation  from  the  potential  temperature  at  700  mb 
and  850  mb  obtained  from  a  nearby  radiosonde  observation.  In  addition,  the  current  algorithm 
uses  A Q  as  input  parameter  for  calculation  of  Tsfc  or  6sfc,  which  can  be  quite  variable  over 
complex  terrain.  The  HRW  model  employs  a  low  level  of  sophistication  and  does  not  invoke  a 
surface  energy  balance  scheme  to  deal  with  surface  temperature.  It  has  also  been  noted  that  the 
model  is  essentially  a  single  surface  layer  model  which  decouples  this  layer  from  the  rest  of  the 
boundary  layer.  This  simplification  imposes  inherent  difficulty  for  the  representation  of  0o. 
Lanicci  (1985)  had  suggested  a  modified  technique  for  calculation  of  0o,  Tsfc,  and  dsfc.  His 
technique,  however,  was  not  verified  by  observations.  Preliminary  tests  of  model  sensitivity  to 
atmospheric  stabilities  (changes  in  A Q  or  b)  have  been  conducted  by  one  of  the  authors  (Huynh). 
Those  tests  show  that  the  HRW  simulation  is  more  sensitive  to  the  change  in  surface  heating  or 
surface  temperature  under  stable  conditions  than  it  is  under  unstable  conditions.  Results  of  those 
sensitivity  studies  will  be  reported  elsewhere. 


(2)  The  results  presented  in  the  previous  subsection  indicate  that  the  HRW  model  simulates  the 
wind  direction  field  with  less  error  than  the  wind  speed  field  for  the  MADONA.  A  brief 
discussion  of  the  relative  error  of  the  wind  direction  (ED)  and  of  the  wind  speed  (As)  appears 
useful.  These  relative  errors  and  its  ratio  (Rds)  are  defined  as 


_  t/D 
~  D 


(45) 
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where  dD  and  dS  denote  the  errors  of  wind  direction  (D)  and  wind  speed  (S),  respectively,  and 

S2  =  if  +  v2  ,  D  =  const.  -(/>,(/)  =  arctaniy/u )  ,  (46) 

where  u  and  v  are  defined  as  eastward  and  northward  component  of  the  wind  vector,  respectively 
according  to  meteorological  convention.  Similarly,  the  relative  errors  in  it  and  v,  and  their  ratio 
(Ruv)  are  defined  as 


_du_  _dv  _E}l 

«  ’  ■‘-'v  ’  Ill  • 

u  v 


A  relationship  between  RDS  and  RIIV  can  be  derived.  Equations  48  and  49  can  be  derived  from  the 
definitions  of  equations  45,  46,  and  47 


(48) 


E»=~^d  arctanf-l  =-^(Eu-Ev). 
U  \u )  Uo 


Equation  50  may  be  derived  from  the  two  equations  above 

R  _  i  {Eu~EMv  _  i  Kv-1)^ 
"  D  u2Eu  +  v2Ev  D  ^Ruv  +tan2 


(50) 


The  wind  direction  (D)  in  equation  50  is  defined  in  the  range:  0  <  D  <  2 n.  For  the  MADONA 
simulations,  D  >  nil  =  1 .0472,  see  figure  4.  Therefore,  the  following  relation  is  true 


( Rm.  - 1 )  tan  (j) 
Ruv  +  tan2  (j)  ’ 


(51) 


Note  that  Ruv  >  0  if  du  and  dv  result  only  from  wind  speed  errors.  A  positive  error  in  wind  speed 
(dS  >  0)  causes  positive  Eu  and  Ev  while  a  negative  error  (dS  <  0)  results  in  negative  Eu  and  Ev. 
Consequently,  the  absolute  value  of  Rds  can  be  estimated  from  equation  5 1 . 

1  -  R  I  forltan^l  >  1 

(52 

1-R;M  for|tan^|<l 


The  above  expression  implies  that  \Rds\  is  <1.0  when  0.5  <  Ruv  <  2.0.  Only  if  Ruv  is  either  large 

or  small,  out  of  the  above  range,  can  the  relative  errors  in  wind  direction  estimation  be  greater 
than  the  relative  errors  in  wind  speed  estimation.  This  situation  is  likely  to  occur  when  the  wind 
direction  is  east,  west  (v  =  0),  south,  or  north  ( u  =  0).  For  the  MADONA  simulations,  such  wind 
directions  occur  infrequently.  Consequently,  the  statistical  results  for  wind  direction  simulations 
appear  better  than  for  wind  speed  simulations.  The  above  argument,  however,  is  not  intended  to 
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provide  a  physical  explanation  for  our  results.  Cauchey  et  al.  (1979),  for  example,  have 
formulated  a  wind  direction  budget  equation.  A  similar  wind  speed  budget  equation  is  not 
difficult  to  derive.  It  is,  however,  beyond  the  scope  of  the  present  study  to  analyze  these  physical 
budgets. 

(3)  Finally,  the  performance  of  the  model  simulations  for  the  strong  wind  speed  category  seems 
undesirable.  There  are  eight  cases  in  this  category,  equation  43c.  They  are  cases  2,  4,  5,  6,  28, 
35,  36,  and  37  (table  2).  As  shown  in  figure  6  (the  lowest  right  panel),  the  datapoints  are 
scattered  around  the  regression  line.  The  correlation  coefficient  (r)  for  the  wind  speed  is  0.1507 
(table  4).  The  difference  between  the  simulated  and  observed  wind  can  be  as  large  as  27  degrees 
in  wind  direction  and  3.5  ms  1  in  wind  speed,  the  extremes  at  an  individual  station. 

A  careful  examination  of  each  of  these  eight  cases  indicates  that  the  wind  fields  at  the  time  of  the 
simulations  all  exhibit  significant  unsteadiness.  There  were  considerable  variations 
(accelerations)  of  wind  vector,  especially  wind  speed,  around  the  time  of  the  simulations,  i.e., 
considerable  8V  /  dt .  If  there  is  a  gusty  wind  under  a  strong  wind  speed  situation,  a  1  ms  1  of 
wind  speed  change  in  a  period  of  5  min,  then  this  fluctuation  alone  is  equivalent  to  an 

—3  —2 

acceleration  on  the  order  of  10  ms  ,  which  can  be  greater  than  the  magnitude  of  acceleration 
induced  by  topography  or  buoyancy.  Recall  that  the  HRW  model  is  not  a  prognostic  model,  just 
a  diagnostic  one,  which  is  applied  to  a  steady  wind  field.  It  comes  with  no  surprise  then  that  the 
model  will  perfonn  less  satisfactorily  for  gusty  (unsteady)  wind  fields. 


4.  Conclusions 


We  have  presented  a  comprehensive  description  of  the  high  resolution  wind  model  which  has 
been  an  undertaking  at  the  U.S.  Army  Research  Laboratory  (formerly  U.S.  Army  Atmospheric 
Sciences  Laboratory)  for  three  decades.  From  the  early  effort  of  model  development,  several 
requirements  have  been  kept  in  mind:  (1)  the  model  can  be  applied  over  very  complex  terrain 
with  limited  meteorological  information  (inputs);  (2)  the  model  is  adaptable  to  small  horizontal 
grid  spacing  on  the  order  of  100  m  for  a  horizontal  domain  smaller  than  10  by  10  km;  (3)  the 
model  is  to  be  run  as  a  combined  system  with  existing  atmospheric  dispersion  and  transport 
models;  (4)  the  model  can  be  used  on  a  small  computer.  Our  HRW  model  uses  Gauss’s  Principle 
of  least  constraint  for  a  variational  adjustment  of  an  initial  estimated  wind  field  in  a  single 
surface  layer  to  conform  with  terrain  structures,  mass  conservation,  and  buoyancy  forces.  This 
basic  approach  and  its  features  appear  to  meet  the  above  requirements. 

The  measurements  from  the  field  study  of  Meteorology  and  Diffusion  Over  Non-Uniform  Areas 
(MADONA)  have  provided  very  valuable  data  to  test  our  HRW  model.  As  shown  in  the  last 
section,  the  results  from  the  present  study  demonstrate  that  the  HRW  model  generally  perfonns 
well  for  a  total  of  39  cases  when  compared  to  MADONA  data.  The  results  of  our  evaluation,  on 
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the  other  hand,  also  indicate  certain  limitations  in  the  HRW  model  applications.  For  example, 
the  model  may  not  simulate  a  gusty  (unsteady)  wind  field  as  successfully  as  a  steady  wind  field. 
Examination  of  these  limitations  and  further  testing  of  the  model  will  allow  us  to  gain  insights 
consistent  with  the  model  physics  and  to  develop  a  next  generation  model. 

Finally,  the  results  of  a  companion  evaluation  of  the  HRW  model  (Williamson  et  al.  2005) 
indicate  that  the  correlation  coefficient  by  itself  is  only  a  weak  indicator  of  the  model  validity.  In 
fact,  that  study  showed  that  the  HRW  model  performance  for  the  MADONA  field  data  is 
actually  slightly  worse  than  the  homogeneous  wind  field  used  to  initiate  it,  and  this  holds  true  for 
all  statistical  measures  that  were  considered. 
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